\defmodule{SmoothingCubicSpline}


Represents a cubic spline with nodes at $(x_i, y_i)$ computed with
the smoothing cubic spline algorithm of Schoenberg \cite{mDEB78a,mPOL93a}.
 A smoothing cubic spline
is made of $n+1$ cubic polynomials. The $i$th polynomial of such a spline,
for $i=1,\ldots,n-1$, is defined as $S_i(x)$ while the complete spline is
defined as
 \[
S(x) = S_i(x), \qquad\mbox{ for }x\in[x_{i-1}, x_{i}].
\]
For $x<x_0$ and
$x>x_{n-1}$, the spline is not precisely defined, but this class performs
extrapolation by using $S_0$ and $S_{n}$ linear polynomials.
The algorithm which calculates the smoothing spline is a
generalization of the algorithm for an interpolating spline.
$S_i$ is linked to $S_{i+1}$ at $x_{i+1}$ and keeps continuity properties for
first and second derivatives at this point, therefore
 $S_i(x_{i+1})=S_{i+1}(x_{i+1})$,
$S'_i(x_{i+1})=S'_{i+1}(x_{i+1})$ and $S''_i(x_{i+1})=S''_{i+1}(x_{i+1})$.

The spline is computed with a smoothing parameter $\rho\in[0, 1]$ which represents its accuracy with respect to the initial $(x_i, y_i)$ nodes. The smoothing spline minimizes
\[
L = \rho\sum_{i=0}^{n-1}{w_i}\left({y_i - S_i(x_i)}\right)^2 +
(1-\rho)\int_{x_0}^{x_{n-1}}\left(S''(x)\right)^2dx
\]
In fact, by setting $\rho = 1$, we obtain the interpolating spline; and
we obtain a linear function by setting $\rho = 0$.
The weights $w_i>0$, which default to 1, can be used to change the contribution
of each point in the error term. A large value $w_i$ will give a large weight
 to the $i$th point, so the spline will pass closer to it.
Here is a small example that uses smoothing splines:

\begin{vcode}

   int n;
   double[] X = new double[n];
   double[] Y = new double[n];
   // here, fill arrays X and Y with n data points (x_i, y_i)
   // The points must be sorted with respect to x_i.

   double rho = 0.1;
   SmoothingCubicSpline fit = new SmoothingCubicSpline(X, Y, rho);

   int m = 40;
   double[] Xp = new double[m+1];       // Xp, Yp are spline points
   double[] Yp = new double[m+1];
   double h = (X[n-1] - X[0]) / m;      // step

   for (int i = 0; i <= m; i++) {
      double z = X[0] + i * h;
      Xp[i] = z;
      Yp[i] = fit.evaluate(z);          // evaluate spline at z
   }

\end{vcode}

\bigskip\hrule
%  \newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{code}
\begin{hide}
/*
 * Class:        SmoothingCubicSpline
 * Description:  smoothing cubic spline algorithm of Schoenberg
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.functionfit;
   import umontreal.iro.lecuyer.functions.*;
   import umontreal.iro.lecuyer.functions.Polynomial;\begin{hide}
import umontreal.iro.lecuyer.functions.MathFunctionWithFirstDerivative;
import umontreal.iro.lecuyer.functions.MathFunctionWithDerivative;
import umontreal.iro.lecuyer.functions.MathFunctionWithIntegral;
import umontreal.iro.lecuyer.functions.SquareMathFunction;
import umontreal.iro.lecuyer.functions.MathFunctionUtil;
\end{hide}

public class SmoothingCubicSpline implements MathFunction,
             MathFunctionWithFirstDerivative, MathFunctionWithDerivative,
             MathFunctionWithIntegral\begin{hide} {

   private Polynomial[] splineVector;
   private double[] x, y, weight;
   private double rho;
\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}

\begin{code}

   public SmoothingCubicSpline (double[] x, double[] y, double[] w,
                                double rho)\begin{hide} {
      if (x.length != y.length)
         throw new IllegalArgumentException ("x.length != y.length");
      if (w != null && x.length != w.length)
         throw new IllegalArgumentException ("x.length != w.length");
      if (rho < 0 || rho > 1)
         throw new IllegalArgumentException ("rho not in [0, 1]");

      splineVector = new Polynomial[x.length+1];
      this.rho = rho;
      this.x = x.clone();
      this.y = y.clone();

      weight = new double[x.length];
      if (w == null) {
         for (int i = 0; i < weight.length; i++)
            weight[i] = 1.0;
      } else {
         for (int i = 0; i < weight.length; i++)
            weight[i] = w[i];
      }

      resolve();
   }\end{hide}
\end{code}
\begin{tabb}
   Constructs a spline with nodes at $(x_i, y_i)$,
   with weights $w_i$ and smoothing factor $\rho$ = \texttt{rho}.
   The $x_i$ \emph{must} be sorted in increasing order.
\end{tabb}
\begin{htmlonly}
   \param{x}{the $x_i$ coordinates.}
   \param{y}{the $y_i$ coordinates.}
   \param{w}{the weight for each point, must be $> 0$.}
   \param{rho}{the smoothing parameter}
   \exception{IllegalArgumentException}{if \texttt{x}, \texttt{y} and
    \texttt{z} do not have the same length, if rho has wrong value, or if
     the spline cannot be calculated.}
\end{htmlonly}
\begin{code}

   public SmoothingCubicSpline (double[] x, double[] y, double rho) \begin{hide} {
      this (x, y, null, rho);
   }\end{hide}
\end{code}
\begin{tabb}
   Constructs a spline with nodes at $(x_i, y_i)$,
   with weights $= 1$ and smoothing factor $\rho$ = \texttt{rho}.
   The $x_i$ \emph{must} be sorted in increasing order.
\end{tabb}
\begin{htmlonly}
   \param{x}{the $x_i$ coordinates.}
   \param{y}{the $y_i$ coordinates.}
   \param{rho}{the smoothing parameter}
   \exception{IllegalArgumentException}{if \texttt{x} and \texttt{y} do
      not have the same length, if rho has wrong value, or if the spline
       cannot be calculated.}
\end{htmlonly}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}

\begin{code}

   public double evaluate (double z)\begin{hide} {
      int i = getFitPolynomialIndex(z);
      if (i == 0)
         return splineVector[i].evaluate(z-x[0]);
      else
         return splineVector[i].evaluate(z-x[i-1]);
   }\end{hide}
\end{code}
\begin{tabb} Evaluates and returns the value of the spline at $z$.
\end{tabb}
\begin{htmlonly}
   \param{z}{argument of the spline.}
   \return{value of spline.}
\end{htmlonly}
\begin{code}

   public double integral (double a, double b)\begin{hide} {
      int iA = getFitPolynomialIndex(a);
      int iB = getFitPolynomialIndex(b);
      double retour;
      int i = 1;

      // Calcule l'integrale
      if (iA == iB) { // les deux valeurs sont sur le meme polynome
         retour = splineVector[iB].integral(a-x[iB], b-x[iB]);
      } else {
         if (iA == 0)
            retour = splineVector[iA].integral(a-x[iA], 0);
         else
            retour = splineVector[iA].integral(a-x[iA], x[iA+1]-x[iA]);
         for (i = iA+1; i<iB; i++) {
            retour += splineVector[i].integral(0, x[i+1]-x[i]);
         }
         retour += splineVector[iB].integral(0, b-x[iB]);
      }
      return retour;
   }\end{hide}
\end{code}
\begin{tabb}Evaluates and returns the value of the integral of the
 spline from $a$ to $b$.
\end{tabb}
\begin{htmlonly}
   \param{a}{lower limit of integral.}
   \param{b}{upper limit of integral.}
   \return{value of integral.}
\end{htmlonly}
\begin{code}

   public double derivative (double z)\begin{hide} {
      int i = getFitPolynomialIndex(z);
      if (i == 0)
         return splineVector[i].derivative(z-x[0]);
      else
         return splineVector[i].derivative(z-x[i-1]);
   }\end{hide}
\end{code}
\begin{tabb} Evaluates and returns the value of the \emph{first} derivative of
 the spline at $z$.
\end{tabb}
\begin{htmlonly}
   \param{z}{argument of the spline.}
   \return{value of first derivative.}
\end{htmlonly}
\begin{code}

   public double derivative (double z, int n)\begin{hide}  {
      int i = getFitPolynomialIndex(z);
      if (i == 0)
         return splineVector[i].derivative(z - x[0], n);
      else
         return splineVector[i].derivative(z - x[i-1], n);
   }\end{hide}
\end{code}
\begin{tabb} Evaluates and returns the value of the \emph{n}-th derivative of
 the spline at $z$.
\end{tabb}
\begin{htmlonly}
   \param{z}{argument of the spline.}
   \param{n}{order of the derivative.}
   \return{value of n-th derivative.}
\end{htmlonly}
\begin{code}

   public double[] getX()\begin{hide} {
      return x.clone ();
   }\end{hide}
\end{code}
\begin{tabb} Returns the $x_i$ coordinates for this spline.
\end{tabb}
\begin{htmlonly}
   \return{the $x_i$ coordinates.}
\end{htmlonly}
\begin{code}

   public double[] getY()\begin{hide} {
      return y.clone ();
   }\end{hide}
\end{code}
\begin{tabb}   Returns the $y_i$ coordinates for this spline.
\end{tabb}
\begin{htmlonly}
   \return{the $y_i$ coordinates.}
\end{htmlonly}
\begin{code}

   public double[] getWeights()\begin{hide} {
     return weight;
   }\end{hide}
\end{code}
\begin{tabb}   Returns the weights of the points.
\end{tabb}
\begin{htmlonly}
   \return{the weights of the points.}
\end{htmlonly}
\begin{code}

   public double getRho()\begin{hide} {
      return rho;
   }\end{hide}
\end{code}
\begin{tabb}   Returns the smoothing factor used to construct the spline.
\end{tabb}
\begin{htmlonly}
   \return{the smoothing factor.}
\end{htmlonly}
\begin{code}

   public Polynomial[] getSplinePolynomials()\begin{hide} {
      return splineVector.clone();
   }\end{hide}
\end{code}
\begin{tabb}
   Returns a table containing all fitting polynomials.
\end{tabb}
\begin{htmlonly}
   \return{Table containing the fitting polynomials.}
\end{htmlonly}
\begin{code}

   public int getFitPolynomialIndex (double x)\begin{hide} {
      // Algorithme de recherche binaire legerement modifie
      int j = this.x.length-1;
      if (x > this.x[j])
         return j+1;
      int tmp = 0;
      int i = 0;

      while (i+1 != j) {
         if (x > this.x[tmp]) {
            i = tmp;
            tmp = i+(j-i)/2;
         } else {
            j = tmp;
            tmp = i+(j-i)/2;
         }
         if (j == 0) // le point est < a x_0, on sort
            i--;
      }
      return i+1;
   }\end{hide}
\end{code}
\begin{tabb}
   Returns the index of $P$, the \class{Polynomial} instance used to evaluate
    $x$, in an \texttt{ArrayList} table instance returned by
  \texttt{getSplinePolynomials()}. This index $k$ gives also the interval in
    table \textbf{X} which contains the value $x$
   (i.e. such that $x_k < x \leq x_{k+1}$).
\end{tabb}
\begin{htmlonly}
   \return{Index of the polynomial check with x in the Polynomial list returned by method{getSplinePolynomials}}
\end{htmlonly}
\begin{code}\begin{hide}

   private void resolve () {
   /*
      taken from D.S.G Pollock's paper, "Smoothing with Cubic Splines",
      Queen Mary, University of London (1993)
      http://www.qmw.ac.uk/~ugte133/PAPERS/SPLINES.PDF
   */

      double[] h = new double[x.length];
      double[] r = new double[x.length];
      double[] u = new double[x.length];
      double[] v = new double[x.length];
      double[] w = new double[x.length];
      double[] q = new double[x.length+1];
      double[] sigma = new double[weight.length];

      for (int i = 0; i < weight.length; i++) {
         if (weight[i] <= 0.0)
            sigma[i] = 1.0e100;
         else
            sigma[i] = 1.0/Math.sqrt(weight[i]);
      }

      double mu;
      int n = x.length-1;

      if (rho <= 0)
         mu = 1.0e100;   // arbitrary large number to avoid 1/0
      else
         mu = 2 * (1 - rho)/(3 * rho);

      h[0] = x[1] - x[0];
      r[0] = 3/h[0];
      for (int i = 1; i < n; i++) {
         h[i] = x[i+1] - x[i];
         r[i] = 3/h[i];
         q[i] = 3 * (y[i+1] - y[i])/h[i] - 3 * (y[i] - y[i-1])/h[i - 1];
      }

      for (int i = 1; i < n; i++) {
         u[i] = r[i-1]*r[i-1] * sigma[i-1]
               + (r[i - 1] + r[i])*(r[i - 1] + r[i]) * sigma[i] +
                 r[i]*r[i] * sigma[i+1];
         u[i] = mu * u[i] + 2 * (x[i+1] - x[i-1]);
         v[i] = -(r[i - 1] + r[i]) * r[i] * sigma[i] - r[i] * (r[i] +
                 r[i+1]) * sigma[i+1];
         v[i] = mu * v[i] + h[i];
         w[i] = mu * r[i] * r[i+1] * sigma[i+1];
      }
      q = Quincunx(u, v, w, q);

      // extrapolation a gauche
      double[] params = new double[4];
      double dd;
      params[0] = y[0] - mu * r[0] * q[1] * sigma[0];
      dd = y[1] - mu * ((-r[0] - r[1]) * q[1] + r[1] * q[2]) * sigma[1];
      params[1] = (dd - params[0])/h[0] - q[1] * h[0]/3;
      splineVector[0] = new Polynomial(params);

      // premier polynome
      params[0] = y[0] - mu * r[0] * q[1] * sigma[0];
      dd = y[1] - mu * ((-r[0] - r[1]) * q[1] + r[1] * q[2]) * sigma[1];
      params[3] = q[1]/(3 * h[0]);
      params[2] = 0;
      params[1] = (dd - params[0])/h[0] - q[1] * h[0]/3;
      splineVector[1] = new Polynomial(params);

      // les polynomes suivants
      int j;
      for (j = 1; j < n; j++) {
         params[3] = (q[j + 1] - q[j])/(3 * h[j]);
         params[2] = q[j];
         params[1] = (q[j] + q[j - 1]) * h[j - 1] + splineVector[j].getCoefficient(1);
         params[0] = r[j - 1] * q[j - 1] + (-r[j-1] - r[j]) * q[j] + r[j] * q[j + 1];
         params[0] = y[j] - mu * params[0] * sigma[j];
         splineVector[j+1] = new Polynomial(params);
      }

      // extrapolation a droite
      j = n;
      params[3] = 0;
      params[2] = 0;
      params[1] = splineVector[j].derivative(x[x.length-1]-x[x.length-2]);
      params[0] = splineVector[j].evaluate(x[x.length-1]-x[x.length-2]);
      splineVector[n+1] = new Polynomial(params);
   }


   private static double[] Quincunx (double[] u, double[] v,
                                     double[] w, double[] q) {
      u[0] = 0;
      v[1] = v[1]/u[1];
      w[1] = w[1]/u[1];
      int j;

      for (j = 2; j < u.length-1; j++) {
         u[j] = u[j] - u[j - 2] * w[j-2]*w[j-2] - u[j - 1] * v[j-1]*v[j-1];
         v[j] = (v[j] - u[j - 1] * v[j-1] * w[j-1])/u[j];
         w[j] = w[j]/u[j];
      }

      // forward substitution
      q[1] = q[1] - v[0] * q[0];
      for (j = 2; j < u.length-1; j++) {
         q[j] = q[j] - v[j - 1] * q[j - 1] - w[j - 2] * q[j - 2];
      }
      for (j = 1; j < u.length-1; j++) {
         q[j] = q[j]/u[j];
      }

      // back substitution
      q[u.length-1] = 0;
      for (j = u.length-3; j > 0; j--) {
         q[j] = q[j] - v[j] * q[j + 1] - w[j] * q[j + 2];
      }
      return q;
   }
}\end{hide}
\end{code}
